{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现基于iris数据集，使用 估计法，建立Sepal.Length、Sepal.Width、Petal.Length对Petal.Width的稳健回归模型"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn import linear_model\n",
    "import random\n",
    "iris=pd.read_csv(\"http://image.cador.cn/data/iris.csv\")\n",
    "x = iris.iloc[:,[0,1,2]].values\n",
    "y = iris.iloc[:,3].values\n",
    "\n",
    "# 对x中的每个变量，随机替换5个较大值\n",
    "row, col = x.shape\n",
    "for k in range(col):\n",
    "    x[random.sample(range(row),5),k] = np.random.uniform(10,30,5)\n",
    "\n",
    "clf = linear_model.LinearRegression()\n",
    "clf.fit(x,y)\n",
    "x_input = np.c_[x,[1]*row]\n",
    "col = col + 1\n",
    "beta = list(clf.coef_) + [clf.intercept_]\n",
    "c = 0.8\n",
    "k = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.041583360442955575\n",
      "0.13288028329048984\n",
      "0.728241189475862\n",
      "0.09522254097415245\n",
      "0.0011120510881178504\n",
      "1.7254532237701597e-07\n",
      "6.460861534568558e-11\n",
      "3.023389358104663e-14\n",
      "6.178605827578578e-17\n",
      "算法收敛，退出循环！\n"
     ]
    }
   ],
   "source": [
    "# 迭代法求解\n",
    "while k < 100:\n",
    "    y0 = np.matmul(x_input,beta)\n",
    "    epsilon = y - y0\n",
    "    delta = np.std(epsilon)\n",
    "    epsilon = epsilon/delta\n",
    "    faik = [0]*col\n",
    "    faikD = np.zeros((col,col))\n",
    "    for i in range(row):\n",
    "        if np.abs(epsilon[i]) <= c:\n",
    "            xi = x_input[i,:]\n",
    "            faik = faik + np.sin(epsilon[i]/c)*xi\n",
    "            faikD = faikD - np.cos(epsilon[i]/c)*np.array([h*xi for h in xi])/(delta*c)\n",
    "            \n",
    "    b = np.matmul(np.linalg.inv(faikD),faik)\n",
    "    beta = beta - b\n",
    "    print(np.max(np.abs(b)))\n",
    "    if np.max(np.abs(b)) < 1e-15:\n",
    "        print(\"算法收敛，退出循环！\")\n",
    "        break\n",
    "        \n",
    "    k = k + 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.4104814942387641"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(np.abs(y - np.matmul(x_input,beta)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.00209922,  0.00210944,  0.41556885, -0.38645105])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现使用普通最小二乘法拟合回归系数，建立线性回归模型"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.09106041273347798,\n",
       " -0.006162626883046202,\n",
       " 0.05853758944345057,\n",
       " 0.3958411864670931]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clf.fit(x,y)\n",
    "list(clf.coef_) + [clf.intercept_]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5670366038984386"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(np.abs(clf.predict(x) - y))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
